---
title: "Further topics in linear regression"
subtitle: "Part 1"
author: "David Barron"
date: "Hilary Term 2018"
fontsize: 10pt
output: 
  beamer_presentation:
    theme: "Madrid"
    toc: false
    slide_level: 3
    keep_tex: true
    df_print: kable
---


```{r setup, echo=FALSE, message=FALSE, warning=FALSE}
knitr::opts_knit$set(width = 100)
knitr::opts_chunk$set(comment =  '')
options('show.signif.stars' = FALSE)
options('digits' = 3)
library(car)
library(tidyverse)
library(haven)
library(effects)
library(ggeffects)
```

# Introduction

### Introduction

This class will meet weekly for 16 weeks.  In odd-numbered weeks I will introduce a substantive topic in quantitative research methods in a 2-hour session.  In even-numbered weeks, I will run "labs", where the focus will be on putting the theory into practice using R.  These sessions will be 1 hour. 

You can find the handouts (including the R code used to generate them) on my personal Weblearn site: https://weblearn.ox.ac.uk/x/MbYn1T.  Datasets that are used in demonstrations and in the labs are also here.  These should be available to anyone with an Oxford single sign on username and password, but let me know if you have any problems.

The class will be assessed by means of an assignment, but that won't be due in until after the end of Trinity Term.  You will be required to carry out an analysis of a dataset of your choice using any of the methods that we have covered in the course.

### Overview of weeks 1-4

- Review of multiple regression
- Modelling
    - Dummy variables
    - Interactions
- Regression diagnostics
    - Normality of residuals
    - Collinearity
    - Model selection
    - Outliers
    - Heteroskedasticity
    - Linearity
- Sample selection bias


### Regression model

$$
y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i}+\dots+\beta_k x_{ki}
+ \epsilon_i,
$$

for $i=1,\dots,n$ sampled observations. 
$\epsilon_i\sim \text{NID}(0,\sigma^2)$.

**Fitted model**


$$
\hat{y}_i = b_0 + b_1 x_{1i} + b_2 x_{2i} + \dots + b_k x_{ki};
$$
$$
 y_i = \hat{y}_i + e_i,
$$
where $b_j$ are estimates of the corresponding $\beta_j$, and the $e_i$
are residuals.  

\emph{Ordinary Least Squares} (OLS) estimates of $b_j$ are those that minimize

$$
\sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n e_i^2
$$

# Modelling

## Dummy Variables

### What are dummy variables?

Often we want to use explanatory variables in regressions that are categorical.  To do this, we have to use _dummy variables_.  Which category a particular observation falls in to is identified by a series of binary (0/1) variables, one fewer variables than there are categories.  That's because there is always one category that does not give us any additional information: if someone isn't a man, they must be a woman and hence we only need a variable identifying whether someone is a man (dummy variable = 1) or isn't a man (dummy variable = 0).  How, though, do we interpret the parameter estimates associated with dummy variables?


### Simulated data

In this example, we have a categorical variable with 4 categories and a continuous variable that are related to a dependent variable in the following way.

$$
y = -.7 x_1 - .2 x_2 + .3 x_3 + .9 x_4 + .4 x_c + \epsilon(0,2)
$$

We first perform a regression of $y$ on the continuous variable, $x_c$, only.


### Regression with continuous variable only

```{r, echo=FALSE, message=FALSE}
dta <- read_csv('dummyvars.csv')

xfac <- ifelse(dta$x1 == 1,'A',
               ifelse(dta$x2 == 1,'B',
                      ifelse(dta$x3 == 1,'C','D')))
dta$xfac <- factor(xfac)
cont <- lm(y ~ xc, data = dta)
summary(cont)
```

### First category excluded

```{r, echo=FALSE}
xf1 <- lm(y ~ xfac + xc, data = dta)
summary(xf1)
```

### Last category excluded

```{r, echo=FALSE}
dta$xfac <- relevel(dta$xfac, 'D')
xf2 <- lm(y ~ xfac + xc, dta)
summary(xf2)
```

### What is the relationship between the two?

Category | $A$ excluded           | $D$ excluded
---------|------------------------|---------------------------
A        | $-0.81$                | $0.92 - 1.73 = -0.81$ 
B        | $-0.81 + 0.70 = -0.11$ | $0.92 - 1.03 = -0.11$ 
C        | $-0.81  + 1.26 = 0.44$ | $0.92 - 0.48 = 0.44$
D        | $-0.81 + 1.73 = .92$   | $0.92$

Parameter estimates give how much that category differs from the _excluded category_.


### Interpreting t-values

Because parameter estimates depend on the arbitrary choice of excluded category, you can't interpret the $t-$values associated with each estimate in the usual way.  To determine whether a dummy variable is statistically signficant, it is conventional to use an $F-$ test, using the formula:

$$
\frac{(RSS_r - RSS_c)/p}{RSS_c/(n-k-p-1)},
$$
where $RSS_r$ is the residual sum of squares  (RSS) from the regression without dummy variables, $RSS_c$ is the SSR from the complete model, $p$ is the number of extra parameters in the complete model, $n$ is the sample size, $k$ is the number of variable in the restricted model (not counting the constant).


### Example

We can get the numbers we want by using the `anova` function in `R`.

```{r}
anova(cont, xf1)
```
You can simply read off the test from this, but you might want to check the formula above using the RSS numbers.

### Duncan's occupational prestige data

This example uses data on occupational prestige. The outcome variable is the percentage of survey respondents who rated an occupation's prestige **excellent** or **good**.  The explanatory variables are _income_, which is the percentage of males in the occupation earning \$3500 or more in 1950; _education_, the percentage of males in the occupation in 1950 who were high school graduates; and _type_, which is a factor distinguishing occupations that are _professions_, _white collar_ or _blue collar_.

### Regression output

```{r, echo=FALSE}
summary(duncan.model <- lm(prestige ~ income + education + type, data = Duncan))
```


### Effect plot

```{r, echo=FALSE}
plot(Effect(c('education','type'), duncan.model), multiline = TRUE, main = 'Effect plot')
```

### Interpretation

You can see that dummy variables have the effect of shifting estimated regression lines up or down.  The lines are parallel to each other.  Here we can see that occupational prestige increases with education and that at all levels of education, estimated prestige is lowest for white collar jobs and highest for professional occupations.

## Interactions

### Motivation

The standard linear regression model implies that the size of the effect of any given explanatory variable on the outcome variable is the same at all values of the other explanatory variables.  What do we do if we think that is not true? For example, the effect of marital status and number of children on wages may be different for men and women.  The standard way of incorporating such interactions is to multiply two variables together:

$$
Y_i = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \beta_3 X_{1i} X_{2i} + \epsilon_i
$$


### Example: Labour Force Survey data

Data from the UK Labour Force Survey gives information about wages as well as age, gender, marital status and number of children.  Wages are tranformed to an hourly basis, and then logged because otherwise they would be very skewed.

```{r, echo=FALSE, message=FALSE, out.width='40%', warning=FALSE}
lfs <- read_dta('C:\\Users\\dbarron\\Dropbox\\Teaching\\MSc teaching\\Advanced Quant\\data/lfs2002.dta')

lfs <- modify_if(lfs, is_character, factor)

#summary(lfs)
lfs$hourpay0 <- ifelse(lfs$hourpay <= 0, NA, lfs$hourpay)
lfs$Loghourpay <- log(lfs$hourpay0)
#lfs$Loghourpay[is.nan(lfs$Loghourpay)] <- NA
#lfs$hourpay[lfs$hourpay <= 0] <- NA

lfs$edage[lfs$edage < 0] <- NA
lfs$allchildren <- with(lfs, numchild + numchil1)

lfs1 <- lm(Loghourpay ~ sex + age + allchildren + married, data = lfs)
lfs2 <- update(lfs1, . ~ . + + sex*married + sex*allchildren)
#summary(lfs1)
ggplot(lfs, aes(x = hourpay0)) + geom_histogram(bins = 35) + xlab('Hourly wage') + theme_bw()

ggplot(lfs, aes(x = Loghourpay)) + geom_histogram(bins = 35) + xlab('Log hourly wage') + theme_bw()
```

### Results without interactions

\footnotesize
```{r, echo=FALSE, size='footnotesize'}
summary(lfs1)
```
All of these estimates are statistically signficant, although the overall model fit is pitiful!


### Results with interactions

\footnotesize
```{r, echo=FALSE, size='scriptsize'}
summary(lfs2)
```


### Interpretation of interactions

```{r, echo=FALSE}
b1 <- coef(lfs1)
b2 <- coef(lfs2)
```

There are two additional parameter estimates, representing the interaction of sex and marital status, and sex and number of children, respectively.  We can work out the effect of being married on log hourly wages for men and women as follows:

Category        |  No interaction           | With interaction
----------------|---------------------------|------------------------
Unmarried women | `r b1[1]`                 | `r b2[1]`
Married women   | `r b1[1] + b1[5]`         | `r b2[1] + b2[5]`
Unmarried men   | `r b1[1] + b1[2]`         | `r b2[1] + b2[2]`
Married men     | `r b1[1] + b1[2] + b1[5]` | `r b2[1] + b2[2] + b2[5] + b2[6]`

You can see that in the first column the difference between being married and unmarried is the same for men and women, but in the second column the differences are much bigger for men than for women.

### Interpretation of interactions 2

```{r, echo=FALSE}
plot(Effect(c('allchildren', 'sex'), lfs2), rug = FALSE, multiline = TRUE, xlab = 'Number of children',
     ylab = 'Log hourly pay', main = 'Effect plot')
```

# Regression diagnostics

## Normality of residuals


### Assumption

The standard assumption of linear regression is that the errors are normally distributed.  If they are not, you will still get unbiased estimates of the regression parameters.  However, the estimates will not (necessarily) be as efficient as they could be (i.e., standard errors will be larger than they need to be).  Hypothesis testing (which relies on us knowing the sampling distribution of estimates) also depends on normality assumption being met.  

### Example

As an example, look at the Labour Force Survey data again but do the regression without taking logs of hourly pay.

\scriptsize
```{r, echo=FALSE, comment='', size='tiny' }
summary(lfs2.nolog <- lm(hourpay0 ~ sex*married + age + sex*allchildren, data = lfs))
```


### Density plot

A density plot is a kind of "continuous histogram".  You can compare the distribution of residuals with a normal distribution with the same mean and standard deviation.

```{r, echo=FALSE, out.width='45%'}
nolog.resid <- residuals(lfs2.nolog)
log.resid <- residuals(lfs2)
n <- length(log.resid)
nl.norm <- rnorm(n, mean(nolog.resid), sd(nolog.resid))
l.norm <- rnorm(n, mean(log.resid), sd(log.resid))
resid.dta <- data.frame(Nolog = nolog.resid, Log = log.resid, Nolog.norm = nl.norm, Log.norm = l.norm)

ggplot(resid.dta) + geom_density(aes(x = Nolog)) + geom_density(aes(x = Nolog.norm), colour = 'red') +
  theme_bw() + ggtitle('Unlogged')

ggplot(resid.dta) + geom_density(aes(x = Log)) + geom_density(aes(x = Log.norm), colour = 'red') +
  theme_bw() + ggtitle('Logged')
```


### QQ-plot

```{r, echo=FALSE, out.width='45%'}
ggplot(resid.dta) + geom_qq(aes(sample = Nolog)) + 
  geom_abline(intercept = mean(resid.dta$Nolog), slope = sd(resid.dta$Nolog), colour = 'red') + 
  theme_bw() + ggtitle('Unlogged')

ggplot(resid.dta) + geom_qq(aes(sample = Log)) +
  geom_abline(intercept = mean(resid.dta$Log), slope = sd(resid.dta$Log), colour = 'red') +
  theme_bw() + ggtitle('Logged')
```

## Multi-collinearity

### Definition

(Multi-)collinearity is the problem of two or more explanatory variables not being independent of each other.  Strictly speaking, this is not a violation of the assumptions of the linear regression model, but when collinearity becomes very high, estimated standard errors become very high and in some circumstances regression parameter estimates can be difficult to obtain.  One way to measure collinearity relies on $R^2_i$, the proportion of the variance of the $i$ th explanatory variable that is associated with the other explanatory variables in the model.  That is, if the regression model is
$$
y = b_0 + b_1 x_1 + b_2 x_x + \dots + b_k x_k + e,
$$
then we regress one explanatory variable on the others:
$$
x_1 = c_0 + c_2 x_2 + \dots + c_k x_k + e
$$
and find the $R^2$ of this second regression.

### Variance inflation factor and tolerance

More commonly, two statistics that are derived from $R^2_i$ are reported.

- **Tolerance** $= (1 - R_i^2);$
- **Variance inflation factor** $= 1 / (1 - R^2_i).$

The VIF is thus the reciprocal of the tolerance.  The VIF (or its square root) is the most commonly reported statistic because it is the impact on the estimated variance (or standard error) of parameter estimates that we are usually most concerned about:
$$
\sigma^2(b_i) = \frac{\sigma^2_\epsilon}{\sum x_i^2} \times \text{VIF}
$$

A common rule of thumb is that a VIF of 10 or above is a source of concern.  However, treat such rules with caution, as it is possible to make matters worse by using common "solutions."

### Example

Calculate the VIF for the Labour Force Survey regression above:
```{r, echo=FALSE}
vif1 <- car::vif(lfs2)
vifm <- data.frame( VIF = vif1)
print(vifm)
#knitr::kable(vifm)
```

You can see that all these VIFs are quite small, so (despite there being two interaction effects, where collinearity can sometimes be a problem), we don't have any concerns about this.  What do we do if there is evidence of high collinearity, though?

### Solutions?

In many cases, there is no straightforward solution; if variables are highly collinear, that's just the way the world is and you can't change it no matter how inconvenient it may be.  For example, it might be difficult to separate the impact of age and years of experience on wages.  It is increasing the risk of failing to reject a null hypothesis even if it is false, so if estimates are significant anyway, you're OK.  If you need to reduce the impact of collinearity, there are a few possibilities.

- Collect more data.  This reduces standard errors, but it may not be practical.
- Combine two or more explanatory variables into a single indicator.  Only an option in (rare) cases where this would make theoretical sense.
- Remove one or more variables from the regression.  This is very risky, and introduces the broader question of how to select the "best" regression model.

### Misspecification bias

```{r, echo=FALSE}
set.seed(1239)
n <- 30
x1 <- rnorm(n, 2, 2)
x2 <- 2 + 0.8 * x1 + rnorm(n, 0, 1)

y <- 4 + 0.5 * x1 - 1.3 * x2 + rnorm(n, 0, 1)
summary(l1 <- lm(y ~ x1 + x2))
```
***

```{r, echo=FALSE}
summary(l2 <- lm(y ~ x2))
```

### Note

This is based on simulated data with $b_0=4$, $b_1=0.5$ and $b_2 = -1.3$. The two explanatory variables are strongly correlated. Removing $x_1$ from the analysis because it is not statistically significant introduces bias in the estimate of $b_2$.


### Things to look out for

- Large change in the parameter estimate of $b_2$ across the two regressions.
- Large change in the $R^2$ across the two regressions.
- **Most important** is your theory; make decisions based on theory, not by blindly following some statistical "rule."
- Consider using one of the step-wise regression methods as an aid to model building.
    - These are particularly appropriate when you are building models with the primary purpose of prediction
    
### Stepwise regression

The basic idea of stepwise regression is to identify a subset of potential explanatory variables that explain as much variance as possible in the outcome variable as parsimoniously as possible.  There are two possible approaches: 

- We start with a minimal model and add variables until there is no improvement in fit;
- We start with all possible variables and remove them until there is no deterioration in fit.

The criterion most commonly used to assess fit is the Akaike information criterion (AIC), which is smaller the better fitting the model, taking account of the number of parameters being estimated.

### Example

Data on credit histories of 1,319 applicants for credit cards.  The outcome variable is the number of major negative reports. _Age_ in years; _Income_ in US dollars/10,000; _Share_ is ratio of monthly credit card expenditure to yearly income; _Owner_ is a factor, whether a home owner; _Dependents_ is number of dependents; _Months_ at current address.

\tiny
```{r, echo=FALSE}
data(CreditCard, package = 'AER')
cc.full <- lm(reports ~ age + income + share + owner + dependents + months, data = CreditCard)

summary(cc.full)
```

### Backwards elimination

```{r, echo=FALSE}
cc.b <- step(cc.full, direction = 'back', trace = 0)
summary(cc.b)
```

### Forwards addition

\scriptsize
```{r, echo=FALSE}
cc1 <- lm(reports ~ 1, data = CreditCard)
cc.f <- step(cc1, scope = ~ age + income + share + owner + dependents + months, trace = 0, direction = 'forward')
summary(cc.f)
```
\normalsize
In this case, both methods give the same answer, which adds to our confidence.  This isn't always the case.  Care needs to be taken using stepwise methods; there is no substitute for thinking!

